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Abstract. Three-dimensional hydrodynamical simula- 
tions are presented for the direct head-on or off-center 
collision of two neutron stars, employing a basically New- 
tonian PPM code but including the emission of gravi- 
tational waves and their back-reaction on the hydrody- 
namical flow. A physical nuclear equation of state is used 
that allows us to follow the thermodynamical evolution of 
the stellar matter and to compute the emission of neu- 
trinos. Predicted gravitational wave signals, luminosities 
and waveforms, are presented. The models are evaluated 
for their implications for gamma-ray burst scenarios. We 
find an extremely luminous outburst of neutrinos with a 
peak luminosity of more than 4 • 10 erg/s for several mil- 
liseconds. This leads to an efficiency of about 1% for the 
annihilation of neutrinos with antineutrinos, correspond- 
ing to an average energy deposition rate of more than 
10 52 erg/s and a total energy of about 10 50 erg deposited 
in electron-positron pairs around the collision site within 
10 ms. Although these numbers seem very favorable for 
gamma-ray burst scenarios, the pollution of the e pair- 
plasma cloud with nearly 10 _1 Mq of dynamically ejected 
baryons is 5 orders of magnitude too large. Therefore the 
formation of a relativistically expanding fireball that leads 
to a gamma-ray burst powered by neutrino emission from 
colliding neutron stars is definitely ruled out. 

Key words: gamma rays: bursts - gravitation: waves - 
elementary particles: neutrinos - stars: neutron - binaries: 
close - hydrodynamics 



1. Introduction 

The merging of two neutron stars that make up a binary 
system is of interest both because it is a powerful source 
of gravitational waves and because it might be the central 
engine of gamma-ray bursts. These mergings have been 
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studied intensely during the past few years, since the oc- 
currence rate is high enough (e.g. Narayan et al. 1991, 
Phinney 1991, Tutukov et al. 1992, Tutukov & Yungelson 
1993, Lipunov et al. 1995, van den Heuvel & Lorimer 1996, 
Lipunov et al. 1997, Prokhorov et al. 1997, Bethe & Brown 
1998) that one might be able to measure the observable 
consequences. However, direct collisions of two neutron 
stars are much rarer [the rate is -C 10~ 10 per year per 
galaxy (Centrella & McMillan 1993) as compared to the 
binary merger rate of 10~ 6 to 10~ 5 per year per galaxy] 
and thus have not received much attention. 

Two previously published sets of hydrodynamical 
models of colliding neutron stars were computed by Cen- 
trella & McMillan (1993) and by Rasio & Shapiro (1992). 
Both used an SPH code to simulate the dynamics of two 
polytropic neutron stars of equal masses with adiabatic in- 
dex 7 = 2 and started at an initial distance of about four 
neutron star radii on parabolic orbits. Using a relatively 
small number of particles (2048), Centrella & McMillan 
(1993) were able to produce a catalog of gravitational wave 
luminosities and gravitational waveforms for different im- 
pact parameters. Rasio & Shapiro (1992) restricted them- 
selves to fewer simulations with higher resolution (16000 
particles) of neutron star coalescence as well as head-on 
collisions. They found that up to 5% of the total mass can 
escape from the systems and that strong shocks occur. 
Their results for the gravitational waveforms and lumi- 
nosities reveal a smaller first maximum (caused by the 
initial free fall of the two stars), followed by the main 
peak associated with the rapid deceleration of the collid- 
ing matter during the propagation of the recoil shocks. 

Neutron star collisions have repeatedly been suggested 
in the literature as possible sources of gamma-ray bursts 
(e.g., Katz & Canel 1996, Dokuchaevfe Eroshenko 1996; 
also Dokuchaev et al. 1998), powered either by neutrino- 
antineutrino annihilation which produces an pair-pho- 
ton fireball, or by highly relativistic shocks which are 
formed during the collision (or coalescence) and eject 
matter at relativistic velocities (Shaviv & Dar 1995). 
Katz & Canel (1996), in particular, developed the idea 
that long gamma-ray bursts might be explained by ac- 
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cretion induced collapse of bare degenerate white dwarfs, 
while short bursts might originate from neutron star col- 
lisions. They estimated the post-collision temperature to 
be kT w 100 MeV, and expected a neutrino pulse with a 
width of 2.5 ms, and a total number of escaping neutrinos 
of 10 58 (refering to Dar et al. 1992). With a mean neu- 
trino energy of about 6 MeV this corresponds to a total 
energy of 10 53 erg emitted in neutrinos. They sketched 
the picture that shortly after the collision, dense lumps of 
hot matter expand and become optically thin to neutrinos 
which are released in a powerful outburst that is luminous 
enough to produce the desired gamma-ray burst energy by 
neutrino-antineutrino annihilation if one assumes an effi- 
ciency of 1%. In order to obtain the observed gamma-ray 
burst rates they had to invoke a total of 10 9 hypothetical 
clusters within z f=a 1 with 10 8 neutron stars each. They 
stated that "we cannot exclude the possibility that such 
clusters are commonly found at the centers of galaxies" . 

Neutron star collisions might also be considered as in- 
teresting, because they could be viewed at as more violent 
variants of the merging process with a more extreme dy- 
namical evolution. This might support the formation of 
relativistic shocks and possibly lead to a larger burst of 
gravitational waves and neutrinos. More violence of the 
merging process can be expected from general relativistic 
effects which were so far neglected in the large majority 
of simulations (see, however, Wilson et al. 1996 for an 
interesting step towards fully general relativistic model- 
ing). Oohara & Nakamura (1997), using a grid-based TVD 
scheme, performed preliminary simulations to compare 
the coalescence with a Newtonian potential to the case 
with post-Newtonian physics. Their neutron stars were 
modeled as two polytropes with 7 = 2, M = 0.62 M©, 
and R = 15 km. The principal difference between the two 
cases was that the post-Newtonian merging indeed turned 
out to be more violent: the impact is more central and 
shocks develop, because "general relativity effectively in- 
creases the gravitational force" . Although the strong shock 
itself has little immediate effect on the gravitational wave- 
form, the dynamics can be changed because of the higher 
temperatures and densities. Wex (1995) and Ogawaguchi 
& Kojima (1996) showed that both spin-orbit and spin- 
spin interactions appear formally to be first order post- 
Newtonian corrections, just as gravitational potential cor- 
rections are (gravitational waves are 2.5 PN), but the in- 
ferred magnitude of these corrections for known compact 
binary systems is actually smaller. 

Our project of simulating neutron star collisions with 
a Newtonian PPM code was motivated by the aspects de- 
scribed in the last two paragraphs. On the one hand, we 
intended to put potential gamma-ray burst scenarios to 
a test, on the other hand we wanted to study a situation 
that mimics the merging of two neutron stars with extreme 
violence and maximal parameters like pre-merging kinetic 
energy and angular momentum in the system. Thus we 
hoped not only to obtain an upper bound on the gravita- 



tional wave emission to be expected from the merging of 
two neutron stars, even with general relativistic effects in- 
cluded. We also wanted to see whether the most extreme 
conditions during the collision of the two stars lead to 
sufficiently large neutrino emission to explain the gamma- 
ray burst energetics by the annihilation of neutrinos and 
antineutrinos emitted during the dynamical event. The 
latter seems impossible in case of the final stages of the 
coalescence of binary neutron stars because the prompt 
neutrino burst, although very luminous, fails by several 
orders of magnitude to produce about 10 51 erg of gamma- 
rays within the short time of only a few milliseconds that 
it takes the two neutron stars to merge into one massive 
body that is most likely going to collapse into a black hole 
on a dynamical timescale (Ruffert et al. 1996). 

The paper is organized as follows. In Sect. the basic 
aspects and new elements of our computational method 
(in extension of Ruffert et al. 1996) and the chosen initial 
conditions for our simulations are given. Our results are 
presented in the following sections, where we describe the 
hydrodynamical and thermodynamical evolution (Sect. ^|), 
the gravitational wave emission (Sect. |j), and the neutrino 
production (Sect. ||) in the colliding stars together with 
the evaluation of our models for the efficiency of neutrino- 
antineutrino annihilation (Sect. |^). Section |?] contains a 
summary and conclusions. 

2. Computational procedure and initial conditions 

In this section we summarize the numerical methods and 
the treatment of the input physics used for the presented 
simulations. In addition, we specify the initial conditions 
by which our different models are distinguished. More de- 
tailed information about the employed numerical proce- 
dures can be found in RJS (Ruffert ct al. 1996) and RJTS 
(Ruffert et al. 1997). 

2.1. Hydrodynamical code 

The hydrodynamical simulations were done with a code 
based on the Piecewise Parabolic Method (PPM) devel- 
oped by Colclla & Woodward (1984). The code is basi- 
cally Newtonian, but contains the terms necessary to de- 
scribe gravitational-wave emission and the corresponding 
back-reaction on the hydrodynamical flow (Blanchet et 
al. 1990). The modifications that follow from the grav- 
itational potential are implemented as source terms in 
the PPM algorithm. The necessary spatial derivatives are 
evaluated as standard centered differences on the grid. 

In order to describe the thermodynamics of the neu- 
tron star matter, we use the equation of state (EOS) of 
Lattimer & Swesty (1991) for a compressibility modulus 
of bulk nuclear matter of K = 180 MeV in tabular form. 
Energy loss and changes of the electron abundance due to 
the emission of neutrinos and antineutrinos are taken into 
account by an elaborate "neutrino leakage scheme". The 
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Table 1. Characterizing parameters and some computed quantities for all models. TV is the number of grid zones 
per dimension, L the size of the largest grid, I the size of the smallest zone, M p< u the gas mass with density 
below 10 11 g/cm 3 at the end of the simulation, and M<j the mass with specific angular momentum larger than 
j* = WKcpior(3i?s)3i? s where R s is the event horizon of a <~ 3M black hole probably forming from the colliding 
neutron stars. M g denotes the mass swept off the grid, M u is the mass leaving the grid unbound, and T ex the maximum 
temperature (in energy units) reached on the grid during the simulation of a model. 
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Table 2. Gravitational- wave and neutrino emission properties for all models. C is the maximum gravitational- wave 
luminosity, £ the total energy emitted in gravitational waves, rh the maximum amplitude of the gravitational waves 
as observed from a distance r, L„ c the stationary value of the electron neutrino luminosity which is reached at 
about 6-10 ms after the start of the simulations, L Pe the corresponding electron antincutrino luminosity, and L Vx the 
luminosity of each individual species of v x (= v^, 9^, v T , D T ). L v represents the total neutrino luminosity at the end 
of the simulation, (e Ve ), (ep e ) and (e Vx ) are the mean energies of the different neutrino and antineutrino flavors. E v p 
denotes the integral rate of energy deposition by neutrino-antineutrino annihilation, averaged over the simulation time 
of 10 ms. 
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energy source terms contain the production of all types 
of neutrino pairs by thermal processes and additionally of 
electron neutrinos and antincutrinos by lepton captures 
onto baryons. The latter reactions act as sources or sinks 
of lepton number, too, and are included as source terms in 
a continuity equation for the electron lepton number. Mat- 
ter is rendered optically thick to neutrinos due to the main 
opacity producing reactions which are neutrino-nucleon 
scattering and absorption of electron-type neutrinos onto 
nucleons. 

More detailed information about the employed numer- 
ical procedures can be found in RJS, in particular about 
the implementation of the gravitational-wave radiation 
and back-reaction terms and the treatment of the neutrino 
lepton number and energy loss terms in the hydrodynam- 
ical code. 

We have extended and improved the numerical treat- 
ment as compared to RJS in several aspects (a comparison 
of published results for coalescing neutron stars obtained 
with the old code against results from the improved one 
will be given in a separate, forthcoming paper): 



(a) Numerical resolution: The presented simulations 
were done on multiply nested and refined grids. With an 
only modest increase in CPU time, the nested grids al- 
low one to simulate a substantially larger computational 
volume while at the same time they permit a higher lo- 
cal spatial resolution of the merged object. The former is 
important to follow the fate of matter that is flung out to 
distances far away from the collision site either to become 
unbound or to eventually fall back. The latter is necessary 
to adequately resolve the strong shock fronts and steep 
discontinuities of the plasma flow that develop during the 
collision. The procedures used here are based on the al- 
gorithms that can be found in Berger & Colella (1989), 
Berger (1987) and Berger & Oliger (1984). Our version is 
described in detail in Section 4 of Ruffert (1992) so only 
the most important features are to be summarized here. 
The individual grids are equidistant and Cartesian with 
each finer grid having a factor of two smaller zone size 
and extent than the next coarser one. Hereby the num- 
ber of zones remains the same for all grids, typically 32 3 
for low-resolution test calculations, 64 3 for our "standard" 
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x axis [ km ] x axis [ km ] 

Fig. 1. Contour plots of Model Ti. (left panels) and H (right panels) showing cuts in a plane containing the x-axis 
which is the symmetry axis of the initial model. The displayed physical quantities are the density together with 
the velocity field (panels a and b) and the temperature (panels c and d). The density is measured in gem -3 , the 
temperature in MeV. The density contours are spaced logarithmically with intervals of 0.5 dex, while the temperature 
contours are linearly spaced, starting with 1 MeV, 3 MeV, 5 MeV, and then continuing with an increment of 5 MeV. 
The bold contours are labeled with their corresponding values (10 10 , 10 12 , and 10 14 gcm~ 3 , and 10, 30, and 50 MeV, 
respectively). In the box in the upper right corner of each panel, the velocity vectors and the time elapsed since the 
beginning of the simulation are given. The mirror symmetry relative to the plane x = and the rotational symmetry 
around the x-axis are broken during the evolution (panels b and d) because of an instability of the contact layer of 
the two neutron stars against shear motions by which numerical fluctuations (panel c) are amplified. 



simulations, and 128 3 for models of special interest where 
high resolution seems desirable. 

(b) Physics input: The table for the Lattimer & Swesty 
(1991) equation of state was extended to higher and lower 
temperatures and now spans 0.01 MeV < T < 100 MeV, 
and also the lower density bound was moved down to now 



Pmin = 5T0 g cm 3 so that the density range now covered 
by the table is 5 • 10 7 gc m - 3 < p < 2.9 • 10 15 gcm- 3 . 
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x axis [ km ] x axis [ km ] 

Fig. 2. Contour plots of Model Ti. (left panels) and H (right panels) showing cuts in a plane containing the x-axis 
which is the symmetry axis of the initial model. The displayed physical quantities are the electron fraction Y e (panels 
a and b) and the entropy (panels c and d), the latter measured in units of Boltzmann's constant k per nucleon. The 
contours of the electron fraction are linearly spaced with intervals of 0.02 below 0.1 and with intervals of 0.05 above, 
the entropy contours are given in steps of lfc/nucleon from 1 to 6, in steps of 2fc/nucleon between 6 and 16, and 
then for values of 20, 25, 30, and 40fc/nucleon. The bold contours are labeled with their corresponding values (0.02, 
0.06, 0.10, 0.20, 0.30, and 0.40 for Y e , and 6, 10, and 20 for the entropy). Maximum values of Y e are above 0.4, of the 
entropy near 30 fc/nuclcon. In the box in the upper right corner of each panel, the time elapsed since the beginning of 
the simulation is given. 



2.2. Evaluation of neutrino- antineutrino annihilation 

In a post-processing step, performed after the hydrody- 
namical evolution had been calculated, we evaluated our 
models for neutrino-antincutrino {yv) annihilation in the 
surroundings of the collided stars in order to construct a 
map showing the local energy deposition rates per unit 
volume. Spatial integration finally yields the total rate 



of energy deposition outside the neutrino emitting high- 
density regions. The "brute force" approach, however, 
which was applied by RJTS, is not feasible any longer be- 
cause it involved explicit summation of the contributions 
of neutrino and antineutrino loss terms of all grid cells and 
at every location where the local annihilation rate was to 
be determined. The computational load for this procedure 
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x axis [ km ] x axis [ km ] 

Fig. 3. Contour plots of Model O showing cuts in the orbital plane for the density together with the velocity field 
(panels a and b) and for the temperature (panels c and d). The density is measured in gem -3 , the temperature in MeV. 
The density contours are spaced logarithmically with intervals of 0.5 dex, while the temperature contours are linearly 
spaced, starting with 1 MeV, 3 MeV, 5 MeV, and then continuing with an increment of 5 MeV. The bold contours are 
labeled with their corresponding values (10 10 , 10 12 , and 10 14 gem -3 , and 10, 30, and 50 MeV, respectively). In the box 
in the upper right corner of each panel, the velocity vectors and the time elapsed since the beginning of the simulation 
are given. 



increases roughly with the third power of the number of 
grid zones if the annihilation map has about the same 
spatial resolution as the grid for the hydrodynamical sim- 
ulation. With the larger number of zones on several levels 
of the nested grid, such a strategy is currently computa- 
tionally impossible. 

Therefore we resort to a different approach which in- 
volves five distinct steps. (1) First, the relevant physical 
quantities are mapped from only that fractional volume of 



the nested grids where most of the neutrino emission (and 
neutrino opacity) comes from, to an equidistant Carte- 
sian grid of fairly high resolution (typically 140 3 ). This 
mapping is done by tri-linear interpolation. (2) Second, 
on this grid, the neutrinosphere for each flavor of neu- 
trino or antineutrino Vi is determined. We define this two- 
dimensional hypcrsurfacc for Vi by the set of those triples 
(x, y, z) where, for each chosen pair of coordinates x and 
y, the vertical optical depth T Vi (z) satisfies the condition 
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Fig. 4. Contour plots of Model O showing cuts in the orbital plane for the electron fraction Y e (panels a and b) and the 
entropy (panels c and d), the latter quantity measured in units of Boltzmann's constant k per nucleon. The contours 
are linearly spaced with intervals of 0.02 for the electron fraction and 1 fe/nucleon for the entropy. The bold contours 
are labeled with their corresponding values (0.02, 0.06, 0.10, 0.20 for Y e and 5 and 10 for the entropy). Maximum 
values of Y e are around 0.16, of the entropy about 10 fc/nucleon. In the box in the upper right corner of each panel, 
the time elapsed since the beginning of the simulation is given. 



T Vi (z) — Az ■ Y^jLz K iti) = 1 where Az is the cell size of 
the Cartesian grid and the local opacity of neutrino 
vi at position j. (3) Thirdly, the local neutrino number 
and energy loss terms of neutrino v t are added up along 
the ^-direction (for each fixed pair (x,y)), and the total 
neutrino emissivity and the corresponding average energy 
of the emitted neutrinos are projected to originate from 
the neutrinosphere of neutrino Vi determined in step (2). 
(4) Fourthly, and most importantly, the energy deposi- 
tion rates by vv annihilation are calculated by integrating 



(summing) only over the two-dimensional neutrinospheres 
as neutrino and antineutrino sources instead of the three- 
dimensional stellar volume as done in RJTS. Additional 
conditions imposed on the construction of the annihilation 
map by RJTS are also used here, i.e., the neutrino emis- 
sion is assumed to occur isotropically around the outward 
pointing local density gradients at the neutrinospheres, 
and the energy deposition by vv annihilation is evaluated 
only in those regions where the baryon density is below 
a certain threshold, typically p < 10 n gcm -3 , and on a 
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Fig. 5. Cut through Model O showing the density con- 
tours together with the velocity field in the plane y = 
perpendicular to the orbital plane near the end of the sim- 
ulation (t — 9.27 ms). The contours are spaced logarith- 
mically with intervals of 0.5 dex, the bold contours corre- 
spond to density values of 10 10 , 10 12 , and 10 14 gcm~ 3 , 
respectively. The velocity vectors are normalized as in 
Fig. |b. 

cylindrical grid with coarser resolution than the Carte- 
sian grid used to represent the neutrino sources, in order 
to limit the costs of the numerically intense calculations. 
(5) Finally, the local energy deposition rates per unit vol- 
ume, E ann (w, <fi, z), are averaged over the </> direction of 
the cylindrical grid: 



E anu (m, z) 



1 

2^ 



E ann (m,(f),z)d(f). 



(1) 



From these average values two-dimensional maps like the 
one shown in Sect. [6] are plotted and integral numbers 
can be obtained by summation along radial or vertical 
directions in the cylindrical grid. 

The total energy deposition rate ("annihilation lumi- 
nosity" ) Lann is obtained from the local values of the en- 
ergy deposition rate per unit volume, £ ann (ro, z), by inte- 
gration over the whole space outside the neutrino emitting 
stellar source: 



E ann (w, (/), z) dV 



(2) 



with dV — m dzu d(f) dz . Given the time dependent func- 
tion L ann (t) one can then calculate the cumulative energy 
deposition by neutrino-antineutrino annihilation accord- 
ing to 



Because the computation is so expensive, however, L ann (t) 
cannot be evaluated on a fine temporal grid, but only at 
a few discrete points in time, ti. With the values L ann (ti) 
we therefore make the following approximation for E ann : 



E, 



1 - 

-Y 



where T{t) is defined by 



T{t) dt 



(4) 



ee ^-{ Ll/B (t)L De (t) [<£„.(<)) + (e*.(t))] + 

2-L^(t)L^(t) [<£„„(*)> + <e P „(t)>]}. (5) 

Here the term multiplied by the factor 2 accounts for the 
equal contributions from v^v^ annihilation and v T v r an- 
nihilation. The form of T{t) in Eq. (|) reflects the main 
dependences of the vv annihilation rate: The energy de- 
position rate increases proportional to the product of the 
neutrino and antineutrino luminosities times the sum of 
the mean energies of the annihilating neutrinos and an- 
tincutrinos; in the denominator the characteristic radial 
extent R(t) of the neutrino source comes from the volume 
integral of Eq. (^) when the latter is performed in spher- 
ical coordinates (compare Eqs. (3) and (10) in RJTS and 
references therein). The ratio appearing in Eq. (^) in the 
sum in front of the time integral then contains geomet- 
rical effects which result from the dependence of the vv 
annihilation rate on the angular distributions of neutrinos 
and antineutrinos. From the hydrodynamical models, the 
neutrino luminosities for the individual neutrino flavors, 
L„ e (i), L& e (t) and the corresponding values for muon and 
tau neutrinos and antineutrinos, are available as functions 
of time as well as the average energies of the emitted neu- 
trinos, (e„ e (t)}, (ee e (i)), and (e Ss (t)) for v^, v^ v T , and 
v T , We found that the typical radial size R(t) of the neu- 
trino emitting object during the phase where by far most 
of the vv annihilation happens is not very strongly time- 
dependent because the wobblings and oscillations change 
the shape and size of the collision remnant only on smaller 
scales but not globally. Therefore instead of T(t) from 
Eq. ([|) we use in Eq. (|j) the simpler expression 



T*(t) ee L Ve {t)L De {t) [(e„.(i)) + (e p .(t))] + 
2-L^(t)L^(t) [K(t)) + K(t))] 



(6) 



E„. 



L ann {t) dt . 



(3) 



Computing E ann from Eq. (H) instead of Eq. (H) involves 
the approximation that the term abbreviated by T{t) or 
T* (t) contains the main time dependence of the integral 
in Eq. (|J). Ideally, the ratio i a nn(^)/-? r *(t) would have 
to be constant. Since this is not the case, we decided to 
employ an average value for a small number N of time 
points where the spatial integral of Eq. (||) was evaluated. 
It turned out that the variation of L ann (t) / ' T* (t) during 
the most interesting phase of the evolution is less than a 
factor 2. 
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2.3. Initial conditions 

We started our simulations with two identical Newto- 
nian neutron stars, each having a baryonic mass of about 
1.63 Mq and a radius of 15 km, which were placed at a 
center-to-center distance of 42 km. The distributions of 
density p and electron fraction Y e = n e /nb (with n e being 
the number density of electrons minus that of positrons, 
and nb the baryon number density) were taken from a one- 
dimensional model of a cold, deleptonized neutron star in 
hydrostatic equilibrium and were the same as in RJS. For 
numerical reasons the surroundings of the neutron stars 
cannot be treated as completely evacuated. The density 
of the ambient medium was set to less than 10 8 g/cm 3 , 
more than six orders of magnitude smaller than the central 
densities of the stars. The total mass on the whole grid, 
associated with this finite density is less than 1O~ 3 M . 

The neutron stars were given the free-fall velocity at 
their respective initial positions (x,y,z) — (—21km, 0,0) 
and (x, y, z) = (21 km, 0, 0). The angle between the veloc- 
ity vectors and the vector connecting the stellar centers 
was varied to produce a head-on collision for Models h, H, 
and TC , and an off-center collision for Models o and O. The 
impact parameter of the latter was chosen to be one neu- 
tron star radius. This impact parameter is the minimum 
distance that two point masses reach along their orbits. A 
compilation of all models together with their characteriz- 
ing grid parameters, initial parameter settings, and some 
results of the numerical simulations is given in Tables 
and|. 

In degenerate matter variations of the temperature 
lead only to minor changes of the internal energy and 
pressure (both are dominated by degeneracy effects) or, 
inversely, the temperature is extremely sensitive to small 
variations of the total internal energy. Therefore any small 
fluctuation caused for example by small numerical errors 
in the calculation of the energy density, will be amplified 
and reflected in temperature fluctuations. Subsequently, 
the neutrino emission, which scales with a high power of 
the temperature T, will be very noisy. For this reason we 
did not start our simulations with cold (T = 0) or "cool" 
(T & 10 8 K) neutron stars as suggested by the inves- 
tigations of Kochanek (1992), Bildsten & Cutler (1992), 
and Lai (1994). Instead, we constructed initial temper- 
ature distributions inside the neutron stars by assuming 
thermal energy densities of about 3% of the degeneracy 
energy density for a given density p and electron fraction 
Y e . The corresponding central temperature was around 
7 MeV, the surface temperature less than half an MeV, 
and the average temperature was a few MeV. Because of 
the small contribution of thermal effects to the pressure, 
these temperatures are unimportant for the neutron star 
structure, and the rapid and violent hydrodynamical evo- 
lution ensures that the results are essentially unaffected 
by the assumed finite initial temperatures. 



The simulations were performed on a Cray-YMP 4/64. 
The models with 64 zones needed about 24 MWords of 
main memory and took approximately 160 CPU-hours 
each, models with 32 zones roughly a factor of 10 less. 
Movies were generated for every model. 

3. Hydrodynamical and thermodynamical 
evolution 

3.1. Head-on collision 

Figures [l] and show the density p, temperature T, elec- 
tron fraction Y e , and entropy s at two different moments 
of the collision: At t = 0.23 ms after the start of the simu- 
lation the largest compression is reached with a maximum 
density of 1.1 • 10 15 gcm~ 3 (see Fig. ||), associated with a 
prominent peak of the gravitational wave luminosity (see 
Figs. [l5| and |l7|), and at t — 2.87 ms when the collision 
remnant has performed several cycles of oscillatory mo- 
tions before it begins to settle into a more quiet state, 
and the neutrino emission starts to decrease from its most 
powerful phase (Figs. |2(] and ^2|). For the earlier time, data 
from the high-resolution Model Ji are plotted in Figs. [I] 
and whereas for the later moment only data from the 
64 3 Model H are available. 

Shortly after the surfaces of the neutron stars have 
touched during their head-on collision (Figs. 0a and[l]3), a 
strong shock wave is generated at the common surface by 
the abrupt deceleration of the matter. This shock prop- 
agates back into the as yet practically undisturbed neu- 
tron star matter. The temperatures directly behind this 
shock reach values of more than 45 MeV, while ahead 
of the shock they were around 5 MeV. The entropy in 
the initially cool neutron star (s <C lfc/nucleon) is in- 
creased to values between lfc/nucleon and 2fc/nucleon 
(Fig. ^p). At the same time, matter is being squeezed out 
perpendicularly to the collision axis and expands behind a 
very strong shock (postshock entropies near 30 k/ nucleon) . 
This shock-heated matter emits electron antineutrinos in 
large numbers (cf. Fig. |2^) and quickly develops from 
initially neutron-rich conditions to a much more sym- 
metric nuclear state characterized by an electron fraction 
Y e 0.4 (Fig. 0a). In contrast, in the interior of the collid- 
ing bodies the composition remains essentially unchanged 
because of the long neutrino diffusion timescales in the 
hot neutron star matter. 

At the collision interface a thin "pancake" like layer 
with very high temperatures up to about 70 MeV (Figs. 0c 
and |) occurs. This sheet is dynamically unstable due to 
shear motions. First indications of a break-down of the 
mirror symmetry relative to the y-z-plane can be already 
seen at t = 0.23 ms in Figs. 0c and ^a. Only a short mo- 
ment later, when the merged bodies bounce back and the 
oblate shape changes into a prolate form, this flat pancake- 
like layer folds asymmetrically and breaks up on a millisec- 
ond timescale (Fig. 0d). Within 3 ms the density distribu- 
tion has already smoothed into a nearly spherical shape 
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Fig. 6. The separation of the density maxima of the two 
neutron stars as a function of time for the three Models H, 
O, and o. 
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Fig. 7. Trajectory described by the density maximum of 
one of the neutron stars in a Eulerian frame for the off- 
center collision, Model O. Different line styles denote dif- 
ferent orbits or phases between closest approaches. The 
kinks are numerical and due to the fact that the position 
of the maximum density is given by the integers corre- 
sponding to the indices of the numerical grid cells. 
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Fig. 8. The maximum density on the grid as a function of 
time for the two Models H and O. 



Fig. 9. The maximum temperature on the grid as a func- 
tion of time for the two Models H and O. 



(Fig. |Tjt>) and most of the kinetic energy of the impact has 
been dissipated by shocks into thermal energy or is carried 



away by ejected matter (see Fig. 12). The collision has lead 



to an increase of the entropy to values near 2 fc/nucleon in 
the merged object (Fig. ^Jc), whereas the shock heated gas 
that forms a very extended, nearly spherical cloud around 



the dense central body, has entropies between 6 /c/nucleon 
and 10 fc/nucleon (Fig. ||d). Ejected clumps of matter with 
even higher entropy (s ~ 20 fc/nucleon) can be identified, 
and positron captures onto neutrons and v e production in 
the hot gas (T ~ several MeV) leads to a rapid increase 
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of the electron fraction to values Y e ~ 0.3-0.4 and higher 
in the expanding debris. 

The formation of shock waves at the moment of the im- 
pact in Model H (Fig. |l|) clearly indicates that the head-on 
collision is strongly inelastic and the dissipation of kinetic 
energy happens very efficiently. Therefore the two neutron 
stars are not able to separate again after the first com- 
pression and reexpansion, however, it takes several (4-6) 
violent oscillations until all the kinetic energy is dissipated 
into heat. The reexpansions produce peaks of the separa- 
tion of the density maxima in Fig. [| while the compres- 
sion phases are reflected in a sequence of very large density 
maxima in Fig. || and temperature maxima in Fig. ^[ The 
steady decrease of the density maxima in Fig. || indicates 
that the oscillations come to a rest within about 4 ms. In 
contrast, the maximum temperature increases from one 
compression to the next (Fig. ^) because of the dissipative 
heating of the stellar plasma. The most extreme temper- 
atures that are reached in Model H during the dynamical 
phase, < t 4 ms, are close to 100 MeV. After settling 
into a static state t > 4 ms, the maximum temperatures 
in the collision remnant is around 40-50 MeV. 

3.2. Off-center collision 

The motion of the two neutron stars is rather complicated 
in case of the off-center collision, Model O. The first phase 
of the infall (t ^ 0.1ms) proceeds essentially along point- 
mass binary orbits (Fig. [?]), until the stars start to touch 
and orbital energy and angular momentum are converted 
into neutron star spin and are consumed by the accelera- 
tion of matter which is flung off the neutron star surfaces. 
The corresponding loss of orbital angular momentum and 
kinetic energy leads to a transformation of the initially 
parabolic orbits into elliptic ones. Even more orbital en- 
ergy is transfered to internal energy when the neutron 
stars come into contact and inelastic interaction sets in. 

After the first closest approach, visible in Fig. || as 
minimum distance <i m i n ~ 3 km of the density maxima of 
the two neutron stars at t sa 0.15 ms, the positions of the 
density maxima describe three nearly elliptic orbits be- 
tween moments of closest approach at t w 0.15ms, 2.4 ms, 
3.3 ms, and 3.6 ms (see also Fig. ||). In Fig. [?] these three 
orbits, represented by the x-y-trajectory of the density 
maximum of one of the neutron stars, are discerned by 
different line styles. The diameters of the orbits (and thus 
the maximum separation of the density maxima) become 
successively smaller due to the inelasticity of the contact 
between the stars during their close encounters. The major 
axis of the first ellipse has a length of more than 30 km, 
corresponding to an apastron separation of the density 
maxima on the grid of about 65 km at t « 1.2 ms (Fig. ||), 
much larger than the initial distance of the neutron star 
centers which was 42 km. This means that the neutron 
stars separate again after their first encounter, but fall 
back towards each other again. Even during the second 



apastron at t ~ 2.9 ms, the density maxima of the two 
stars have a distance of about 35 km, which is larger than 
the sum of their initial radii 2i? nSi i = 30 km. During the 
third and the subsequent quasi-elliptic orbits the neutron 
stars are not able to separate again. While the first bound 
orbit has a period of about 2.2 ms, the second, smaller 
one has only ^0.9 ms, and the following are even shorter 
(Figs | and 0). 

In Figs. [I and |] the density distribution, temperature, 
electron fraction, and entropy per nucleon are plotted for 
Model O in the orbital plane at two different stages dur- 
ing the off-center collision. The left panels show the results 
for a time close to the apastron of the first orbital ellipse, 
i.e., a little more than a millisecond after the first clos- 
est approach. The neutron stars are tidally strongly de- 
formed and gas has been swept into the surroundings dur- 
ing the first direct contact and interaction (Fig. ||a). There 
is a dense gas bridge (p ~ 10 12 -10 13 g cm~ 3 ) between the 
stars which continue to wobble and oscillate along their 
orbits. The temperature has climbed to nearly 40 MeV in 
distinct hot spots where the gas bridge hits the denser 
cores of the stars (Figs. ||c and ^|), whereas the extended 
cloud of gas surrounding the orbiting bodies has a tem- 
perature of 1-5 MeV. In this ambient gas, v e production 
by positron capture onto neutrons has raised the electron 
fraction from initially less than about 0.05 to maximum 
values around 0.2 (Fig. ^a). The maximum entropy val- 
ues of s ~ Ilk per nucleon are produced by bow shocks 
in front of the rather dilute (p ~ 10 9 -10 10 g cm -3 ) clouds 
reaching outward from the two neutron stars. Clumps of 
gas with entropy s ~ 5-6 k / nucleon are scattered in the 
surroundings of the collision site. Note that in contrast to 
the head-on collision, there is no shock heating of the inte- 
rior of the neutron stars. Up to the end of the simulation 
the high-density cores of the neutron stars retain their low 
initial entropy, even after they have merged into one body. 

The right panels in Figs. || and |J show snapshots at 
a time near the end of the simulation (t = 9.27 ms). 
The distributions of density, entropy, and electron fraction 
have become roughly circular in the x-y-plane: A compact, 
dense central body (p > 10 12 gem -3 ) with T ~ 5-10 MeV 
(outside of two distinct hot spots where T ~ 40 MeV), 
Y e ~ 0.04-0.1, and s ^ 7fc/nucleon is surrounded by an 
extended envelope with somewhat larger Y e and s (but 
lower temperature), which is rapidly rotating and which 
is stabilized by centrifugal forces. Therefore the vertical 
extension of the gas envelope is significantly smaller than 
its diameter in the orbital plane; the density contour corre- 
sponding to p = 10 10 g cm -3 extends to a radius of about 
130 km in the orbital plane, whereas its butterfly shape 
has a maximum vertical height of roughly 70 km (Fig. |s|). 
Even the compact core is rotationally deformed with an 
axis ratio of 1:1.5. Nevertheless, there is a considerable 
amount of gas at large heights \z\ above and below the or- 
bital plane. By spatial integration we find 3.8 x 10~ 2 Af© 
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Fig. 10. Cumulative amount of matter flowing off the nu- Fig. 11. Cumulative amount of matter that is unbound 
merical grid (measured in 10 -2 Mq) as a function of time when it leaves the numerical grid (in units of Mq, plot- 



from the beginning of the simulations. 
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ted logarithmically) as a function of time for Models o, 
O, h, and H. The unbound mass is determined by the cri- 
terion that the sum of the specific gravitational, kinetic, 
and internal energies is positive. 
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Fig. 12. The internal, kinetic, and potential energies of 
the matter on the grid and the sum of these energies as 
functions of time for Model H. Also plotted is the cumu- 
lative energy carried away in gravitational waves. 



Fig. 13. Same as Fig. |l| but for Model O. The thin solid 
line marks the E = level. 



for \z\ > 40 km, 2.2 x 10~ 2 Mq for \z\ > 50 km, and still 
6.0 x 10~ 3 Mq for \z\ > 80 km. 

3.3. Comparison of head-on and off-center collisions 

The different dynamical evolutions of the head-on colli- 
sion, Model H, and the off-center collision, Model O, are 



reflected in the different time histories of the maximum 
density (Fig. ||) and the maximum temperature (Fig. |9|) 
on the grid. Model H shows very large amplitudes of the 
maximum density at the moments of strongest compres- 
sion with up to a factor of ~ 10 larger values compared 
to the return points of expansion phases. These density 
fluctuations are damped in a sequence of 5-6 strong os- 
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Fig. 14. Total angular momentum (z-component) of the 
matter on the grid (upper curves) and cumulative angu- 
lar momentum loss by gravitational wave emission (lower 
curves) as functions of time for Models O and o. Gravita- 
tional waves carry away about 13% of the initial angular 
momentum of the colliding neutron stars and another non- 
negligible amount (~ 3.5%) is taken away by the matter 
flowing off the numerical grid. 

cillations, and a stationary value is reached after about 
4 ms. The temperature evolution reveals spikes and val- 
leys, respectively, at the same moments, but there is a 
general trend of the temperature to increase during the 
dissipative oscillatory compressions and reexpansions. Al- 
though in Model O the variations of the separation of 
the density maxima are much more pronounced than in 
Model H (Fig. |), the maximum density and temperature 
on the grid show much less extreme fluctuations because 
the neutrons stars describe orbits around each other and 
do not crash violently into each other. One can recognize 
peaks of the maximum density correlated with the mo- 
ments of closest approaches, t w 0.15 ms, 2.4 ms, 3.3 ms, 
and 3.6 ms (compare Fig. || and Fig. [?]). The whole evolu- 
tion of Model O is much less violent than that of Model H. 
Nevertheless, on a longer timescale t <; 4 ms, both mod- 
els settle to roughly the same maximum temperature of 
40-50 MeV. The maximum density in Model O becomes 
even somewhat larger than in Model H, because the latter 
has been heated to higher entropies and therefore thermal 
pressure causes an expansion of the collision remnant. In 
addition, the compact core of the remnant of Model H is 
less massive as a result of the larger mass loss during the 
collision. 

The more violent collision and therefore higher tem- 
peratures in Model H push more matter off the grid than 
in Model O (Fig. FlOf). Also a larger fraction of this matter 
gets unbound (Fig. [Ll]) which is the case when the total 
specific energy of the gas, defined as the sum of the spe- 



cific kinetic, internal, and gravitational potential energies, 
becomes positive. In Model H about 1.5 x 10~ 2 M Q are 
able to escape the gravitational potential of the collision 
remnant, whereas it is little more than one tenth of this 
amount in Model O (cf. Table |l|) . Obviously, the angular 
momentum of Model O (see Fig. [l4]) and the associated 
centrifugal shredding of material can hardly compete with 
the ejection of gas in the strong shock waves occurring in 
Model H. We note in passing that Rasio & Shapiro (1992) 
(RS) found a significantly larger amount of mass loss (up 
to about 5% of the total mass) in their simulations of 
head-on collisions between two identical 7 = 2 polytropes, 
compared to only ~ 0.46% that can escape from the sys- 
tem in our Model H. The difference is presumably caused 
by a combination of reasons, the use of different equations 
of state and correspondingly different stellar structure and 
mass (Lattimer & Swesty nuclear EOS here vs. 7 = 2 adi- 
abatic EOS by RS), the inclusion of gravitational wave 
back-reaction on the hydrodynamical flow in our simula- 
tions, and last but not least the use of different numerical 
schemes (Eulerian PPM here vs. SPH by RS) in combi- 
nation with possibly different criteria to determine the 
unbound mass. 

The temporal evolutions of the internal, kinetic, and 
gravitational potential energies (Figs, [l^ and |l^) show 
structures that correspond to the stages of the dynami- 
cal interactions of the colliding stars. The internal energy 
of Model H oscillates strongly with maxima at the mo- 
ments of strongest compression (compare Figs. || and |l2|) 
which coincide with maxima of the kinetic and minima of 
the potential energy. There is a general trend for the in- 
ternal energy to increase with time. This corresponds to a 
decrease of the kinetic energy (while the potential energy 
fluctuates around an essentially constant level) and thus 
reflects the action of dissipative forces. In the off-center 
collision of Model O the kinetic energy is much less effi- 
ciently converted into internal energy. The latter exhibits 
a continuous increase with superimposed, but much less 
dramatic, local peaks at the instants of closest approach, 
also coinciding with maxima of the kinetic and minima 
of the potential energy. This reflects the dynamical trans- 
formation between the different energy forms. Near the 
end of the computed evolution, practically all of the ki- 
netic energy of Model O is rotational energy, E^in = E Iot • 
Using the values from Fig. [l^, we determine the following 
ratio of kinetic energy to the gravitational binding energy: 
p = E kin /E rot » 0.08. 

The total energies of Model H and Model O (defined 
as the sum of kinetic, internal, and potential energies of 
all gas on the grid) are also shown in Figs. [l2| and O, re- 
spectively. In Model H minor variations of the total energy 
between t « 0.2 ms and t ~ 2.2 ms are numerical because 
of the extremely violent collision. Since the energy carried 
away by gravitational waves and mass loss off the numer- 
ical grid is negligible, the total energy at the beginning 
and end of the simulation are nearly equal. In contrast, in 
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Model O the energy emitted in gravitational waves leads 
to a gradual decrease of the total energy of the gas on the 
grid. 

Fig. |l4| shows the total angular momentum (z-com- 
ponent) of the gas on the grid and the cumulative value 
of the angular momentum carried away by the emitted 
gravitational waves as a function of time for Models O 
and o. There is very good agreement of both calculations 
concerning the angular momentum loss in gravitational 
waves. This indicates that the overall mass distribution 
in the neutron stars (which enters the calculation of the 
mass quadrupole moment needed for the evaluation of the 
gravitational-wave source terms) is sufficiently well rep- 
resented even on the coarser grid of Model o. Most of 
the decrease of the gas angular momentum is explained 
by the gravitational wave emission. An additional effect 
comes from the mass loss off the computational grid at 
4 ms t 8 ms (Fig. |l^) which removes an angular mo- 
mentum of about AJ Z ss M^r g v g rj 2 x 10 48 gcm 2 s _1 
(with M g taken from Table |] and r g ss 160 km being the 
grid radius and v g ~ 3.5 x 10 9 cms" 1 the nearly tangential 
velocity of the gas when it leaves the grid) or 3.5% of the 
initial angular momentum. However, although the gravi- 
tational wave loss and the mass flowing off the grid are 
very similar in both models, Model o exhibits a steeper 
decrease of the total angular momentum at t ^ 3 ms than 
Model O. This difference is purely numerical and caused 
by the coarser grid resolution of Model o. Even in the bet- 
ter resolved calculation, Model O, about 7% of the initial 
angular momentum are destroyed by numerical effects at 
the end of the simulation. 

The relativistic rotation parameter is defined as a = 
Jc/(GM^ ot ) where M to t is the total mass of the system 
(Mtot = 2M initially with M being the mass of one of 
the neutron stars) and J is the total angular momentum 
relative to the center of mass. We have an initial value 
of a = 0.60 and find a value of a > 0.47 (in Model O) 
at the end of the simulation after angular momentum has 
been removed from the system by gravitational waves and 
ejected matter. Since a < 1 rotation seems unable to pre- 
vent the collapse of the collision remnant to a black hole 
if the remnant mass exceeds the maximum stable mass of 
the employed equation of state (see also RJS and Rasio 
& Shapiro, 1992, and references therein). Thermal pres- 
sure can increase this stable mass limit only insignificantly 
(cf. Goussard et al. 1997) and, if so, only during the tran- 
sient period of neutrino cooling (note that the interior 
of the neutron stars retains its initial low entropy in the 
off-center collision, see Figs. ||c and d, and thus the tem- 
perature remains fairly low), and also rotation leads to an 
increase of the upper stability limit on the baryon mass 
by only ^ 20% (Friedman et al. 1986; Friedman & Ipser 
1987). 



4. Gravitational waves 

The gravitational-wave luminosities and the cumulative 
energy loss in gravitational waves as functions of time for 
Models H and h and Models O and o are shown in Figs, [if] 
and [IB], respectively, and the corresponding gravitational 
waveforms h + and h x for Models H and O are plotted in 
Fig. 

In Model H the most prominent luminosity spike is 
created at the moment when the two neutron stars crash 
into each other and the gas flow is abruptly decelerated 
and redirected by the recoil shocks (t « 0.22 ms, cf. Fig. |lj) 
which leads to a rapid change of the mass quadrupole mo- 
ment. The peak luminosity reaches about 3.7 x 10 55 erg/s 
(see also Fig. |l7]) . A precursor with about 1 /4 of the max- 
imum luminosity is caused by the increasing tidal defor- 
mation of the neutron stars as they approach each other. 
After this initial outburst the gravitational- wave luminos- 
ity continues to oscillate regularly with a period between 
two maxima of roughly half a millisecond but with peaks 
at least one order of magnitude below the maximum lumi- 
nosity. This indicates that the bulk of the matter quickly 
adopts a more or less spherical distribution. Within about 
8 ms the luminosity falls by more than 5 orders of mag- 
nitude to less than 10 50 erg/s. This dramatic drop is re- 
flected in the gravitational waveforms of Model H which 
indicate that after ~ 4 ms the activity has essentially 
ceased. This coincides with the complete dissipation of 
the kinetic energy at that time (see Fig. [l2]). More than 
50% of the total energy emitted in gravitational waves are 
contained in the luminosity spike and after about 3 ms 
only insignificant further contributions are added (cf. cu- 
mulative energy loss in Fig. |l5|). From Figs. |l5| and [l7|onc 
learns that the coarser resolution of Model h leads only 
to minor differences compared to Models H and Ti with 
the tendency to overestimate the gravitational-wave lumi- 
nosity and emitted energy. Models H and TL are hardly 
distinguishable in Fig. [l7| and the computations seem to 
be converged. 

Model O is an approximately ten times more ener- 
getic source of gravitational waves than Model H and 
emits a total energy of 2 x 10~ 2 Mqc 2 . On the one hand 
the first luminosity maximum at f sj 0.22 ms is nearly 
twice as high (~ 6 x 10 55 erg/s) as in Model H and nearly 
four times longer (half width about 0.2 ms compared to 
0.05ms) (Fig. 18) On the other hand Model O continues 
its strong emission of gravitational waves for the whole 
computation period of 10 ms during which the luminos- 
ity on average stays around 10 erg/s and hardly ever 
drops below 10 53 erg/s. This can be explained by the rapid 
change of the quadrupole moment of the system as the two 
neutron stars repeatedly come close and separate again 
on their quasi- elliptic orbits around each other (cf. Fig. |7|) 
and by the large kinetic energy retained as rotational en- 
ergy of the collision remnant at the end of the simulation 
(Fig. |l^). Several peaks of the gravitational- wave luminos- 
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Fig. 15. Gravitational-wave luminosity and cumulative Fig. 16. Same as Fig. [15] but for Models O and o. 
energy emitted in gravitational waves (measured in units 
of 10~ 2 Mqc 2 ) as functions of time for Models H and h. 




Fig. 17. Gravitational- wave luminosity as a function of 
time for the head-on collision Models Tt, H and h, com- 
pared with the results from Centrella & McMillan (CM, 
1993, Fig. 8) and from Rasio & Shapiro (RS, 1992, Fig. 9). 



Fig. 18. Gravitational-wave luminosity as a function of 
time for the off-center collision Models O and o, compared 
with the result from Centrella & McMillan (CM, 1993, 
Fig. 11). 



ity can be correlated with the moments of closest approach 
of the density maxima of the two neutron stars (compare 
Figs. |l6] and ^), and two strong, short minima of the lu- 
minosity (at t ~ 1.2 ms and t « 4.1ms) coincide with 
instants of maximal separation. In Model O only about 
30% of the total energy emitted in gravitational waves is 
contained in the first luminosity spike, another ~ 50% are 
added during the second and third periastrons (t rs 2.4 ms 
and t w 3.3 ms), and a non- negligible fraction (~ 20%) 
comes at later times. The waveforms of Model O (Fig. [19|) 
exhibit a rather irregular structure during the first 4-5 ms 



of the evolution with a peak amplitude at t « 0.22 ms 
which is about twice as high as in Model H. They continue 
with a very regular, slowly damped sinusoidal modulation 
until the end of the simulation at t w 10 ms. 



Finally, in Figs. |L7| and nj we compare the first large 
spike of the gravitational-wave luminosity with the corre- 
sponding results obtained by Centrella & McMillan (1993) 
and Rasio & Shapiro (1992). In order to do that we rescale 
and renormalize their dimensionless quantities to physical 
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Fig. 19. Gravitational waveforms, h + and h x , for Mod- 
els H and O, respectively, as functions of time. Note that 
the gravitational- wave field h x with the cross polariza- 
tion is practically zero in case of the head-on collision, 
Model H, because there are only very small deviations 
from the mirror symmetry relative to the y-z plane when 
the hot contact layer of the two neutron stars becomes 
unstable against shear motions (compare Figs. |l| and |J). 



units for our chosen neutron star parameters by using the 
dynamical timescale 



\GM ) 



1/2 



0.125 ms 



(7) 



with M — 1.63 Mq as the mass and R = 15 km as the ra- 
dius of the neutron stars, and by the scaling the luminosity 
with the factor 

5 



U= - 



1 / GM 
G \~^R 



3.86 • ltt>* erg/s 



(8) 



Moreover, there is a time shift between their calculations 
and ours because of the different initial center-to-center 
distances of the neutron stars. While it was ao = 42 km 
in our models, Rasio & Shapiro (1992) and Centrella & 
McMillan (1993) assumed a x = AR = 60 km for their 
head-on collisions, and Centrella & McMillan (1993) took 
ai 4.53i? « 68 km for the off-center case. The time 
lag for the parabolic infall trajectories of the two masses 
with different initial separations can be determined as the 
difference of the times to reach the minimum distance (pe- 
riastron), to be (Roy, 1982, Eqs. (4.82) and (4.85)) 



At = 



3VGM~ t 



(9) 



with Ai ee (a, + 2R p )- s /a i - R p (i = 0, 1). M t = M x + 
Mi = 2M is the total mass of the system and the perias- 
tron distance R p is the separation of two point particles 
of mass M at closest approach, which is set to zero for the 
head-on collision and to R p = 1R for the off-center case. 

Taking into account these aspects in Figs. [l7] and 
we find good overall agreement between the gravitational- 
wave luminosities from the different calculations, both in 
shape and magnitude. The remaining minor discrepancies 
can be attributed to the use of the Lattimer & Swesty 
(1991) EOS in our calculations instead of an adiabatic 
EOS with constant index 7 = 2, the possible influence 
from the inclusion of gravitational- wave back-reactions in 
our models, and the effects resulting from different numer- 
ical schemes and different resolution. 



5. Neutrino emission 

Because of the different dynamical evolution, the head-on 
and off-center collisions also show distinctive differences 
in the neutrino emission. The total neutrino luminosities 
for Models H and h and Models O and o as functions of 
time are plotted in Fig. 



20 



One can see that in the head-on collision (Models H 
and h) a first very luminous burst of neutrinos with a peak 
flux of more than 5 x 10 53 erg/s is produced at the moment 
when the two stars crash into each other and the neutron 
star matter is shock heated and hot gas is squeezed out 
perpendicular to the collision axis (see Fig. [l]a and c). 
A second and third luminosity maximum, however, with 
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Fig. 20. Total neutrino luminosities (sums of the con- 
tributions from all neutrino and antincutrino flavors) as 
functions of time for Models H and h and Models O and 
o. 

neutrino radiation, model H 
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time [ms] 

Fig. 22. Luminosities of the individual neutrino and an- 
tineutrino flavors (y e , D e , and the sum of all v x ) as func- 
tions of time for Model H. 

about 8 times larger peak fluxes of more than 4x 10 54 erg/s 
are present in the time interval 1.5 ms ^ t 3.5 ms which 
is the time of maximum temperatures in the collision rem- 
nant (Fig. [9]) when the hot matter starts to expand and 
to spread out over a larger volume (Figs. [l]b and |l|d). A 
luminosity of L v = 4 x 10 54 erg/s « ^Rl\{i ■ |a rad T 4 ) 
corresponds to a neutrinosphere with radius R v ~ 50 km 
which radiates neutrinos and antineutrinos of all flavors as 
black body with temperature T v w 8.5 MeV. Figures |l|b 
and d confirm that the massive, nearly spherical central 



50 




O F ... 1 ... 1 ... 1 ... 1 ... 1 . = 
2 4 6 8 10 

time [ms] 

Fig. 21. Average energies of emitted neutrinos v ei v ei and 
v x (= v^, Pfj,, v T , v T ) as functions of time for Models H 
and O. 



neutrino radiation, model O 




body of the collision remnant, which is embedded in a 
cloud of less dense and cooler gas, has a radial size and 
"surface" temperature in this estimated range. The dura- 
tion of this extremely luminous burst is only about 2 ms, 
after which the total luminosity settles down to a much 
lower but still sizable value around 10 54 erg/s. During the 
~ 10 ms of simulation time, Models H and h emit an en- 
ergy of 6.7 x 10~ 3 Mqc 2 or 1.2 • 10 52 erg in neutrinos, of 
which 50% are radiated during the double peak of the lu- 
minosity. Model H is a stronger source of neutrinos than 
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of gravitational waves; the latter carry away only about 



emitted radiation 



2.2 x 10" 3 M Q c z (Fig. 



Nevertheless, despite of the ex- 



tremely high luminosity the total energy radiated in neu- 
trinos by Model H is still one order of magnitude below 
the estimates by Katz & Canel (1996). 

During the break-out of the bounce shock at t w 0.2 ms 
neutrinos with average energies in excess of 50 MeV are 
emitted from Model H (Fig. [2l]). A second phase of very 
high mean energies coincides with the two luminosity 
spikes and therefore the phase of highest temperatures in 
Model H. As is well known from type II supcrnovae (see, 
e.g., Janka 1993), neutron-rich, hot neutron star matter is 
more opaque to v e than to v e because of frequent captures 
of the v e on the abundant free neutrons. Heavy-lepton neu- 
trinos (y x = Vp, Up, v T , D T ) are even less strongly coupled 
to the stellar medium since their opacity is dominated by 
neutral-current neutrino-nucleon scatterings but they do 
not interact with nucleons via charged-current reactions. 
For these reasons v e decouple energetically from the hot 
plasma at higher densities and thus usually higher tem- 
peratures than v e , and v x at even higher densities and 
temperatures. This explains why typically the mean en- 
ergy of the emitted v e is lower than that of v e which in 
turn is below the average energy oiv x (see Fig.|l]). After 
t <; 4ms we obtain mean energies of (e„ e ) w 10-13 MeV, 
(e Pe ) fa 15-20 MeV, and (e„J » 20-25 MeV which is in the 
range of values found during the neutrino cooling phase 
of newly formed neutron stars in type II supernovae. De- 
spite of this generic ranking of the mean energies, the v e 
luminosity of the collision remnant is larger than the lumi- 
nosity in each individual type of v x , and the v e luminosity 
of the neutron-rich, hot neutron star matter dominates the 
v e luminosity (Fig. |22| ). The difference in ranking between 
neutrino luminosities and mean energies of emitted neu- 
trinos reflects the fact that the emission is not like an ideal 
black-body, but non-equilibrium effects play a role. More- 
over, Fig. 25 shows that there is an extended region (with a 
broad range of temperatures and densities) where the neu- 
trino fluxes are fed by local neutrino energy losses. There- 
fore the neutrino emission can hardly be characterized by 
the conditions of thermodynamical equilibrium at a well 
defined neutrino emitting surface ( "neutrinosphere" ) as- 
sociated with each type of neutrino or antincutrino. 

In the off-center collision, Models O and o, the neutrino 
luminosity reveals a steady increase and does not have 
such pronounced maxima as seen in Model H (Fig. poj ), 
although some fluctuations up to a factor of 2 are present. 
At the end of our simulations (t ss 10 ms), Model O has 
a total neutrino luminosity of about 7 x 10 53 erg/s which 
is only 30% less than in Model H. Because of the lack 
of a phase of extremely high neutrino emission and the 
strong production of gravitational waves during its whole 
evolution, Model O loses only half the energy in neutri- 
nos as Model H but is a much stronger gravitational-wave 
source (Fig. 24). A comparison of the two plots in Fig. ^5] 
shows that in Model O the neutrino emitting gas cloud 




4 6 
time [ms] 



Fig. 24. Comparison of the cumulative energies emitted 
in neutrinos and gravitational waves (in units of Mqc 2 ) 
as functions of time for the four Models H, h, O, and o. 



(at t w 10 ms) is more extended than in Model H (at 
t ~ 2.5 ms) but the peak values of the local energy loss 
rate in neutrinos are only around 3 x 10 32 ergcm~ 3 s _1 in 
Model O whereas they are above 10 34 erg cm -3 s" 1 in case 
of Model H. The mean energies of the emitted neutrinos 
(Fig. |l|) are very similar in Models O and H, and also the 
relative contributions of the different neutrino types to 
the energy loss are roughly similar, as can be seen from a 
comparison of the relative sizes of the individual neutrino 
luminosities in Figs. E3| and G2. 



6. Neutrino-antineutrino annihilation 

The rate of energy deposition by neutrino-antineutrino an- 
nihilation increases, roughly, with the product of the lo- 
cal neutrino and antineutrino energy densities times the 
mean energy of these neutrinos times a factor that ac- 
counts for the angular distribution of the neutrinos (the 
process is very sensitive to the angle at which neutrinos 
and antineutrinos collide, for details, see RJST). In the 
spherically symmetric situation this can be converted into 
a product of the neutrino and antineutrino luminosities 
times the sum of the mean neutrino and antineutrino en- 
ergies times a normalized factor which results from the 
phase space integration over the local neutrino distribu- 
tion functions and which depends on the geometry of the 
considered problem. This was used to arrive at the approx- 
imate description summarized in Eqs. (P])~@ of Sect. 2.2 
employed here to evaluate Model H for the energy deposi- 
tion by vv annihilation in the surroundings of the collision 
remnant. 

Neutrino-antineutrino annihilation is considered as a 
mechanism to pump energy into a fireball consisting of e + , 
e~ , photons, and a small number of baryons. This fireball 
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Fig. 25. Contour plots of the local energy loss rates due to the emission of neutrinos and antineutrinos of all flavors. 
The left plot shows a cut in the z = plane for Model H at time t — 2.47 ms, the right plot displays a cut in the orbital 
plane of Model O at time t — 9.94 ms. Note that the visualized region is larger in the latter figure. The contours 
are spaced logarithmically in steps of 0.5 dex, the bold lines are labeled with their corresponding values in units 
ergcm _3 s _1 . The dashed lines indicate the approximate positions of the neutrinospheres (defined where the optical 
depths for v e , v e , and v x , respectively, are approximately unity). 



was suggested to be a possible source of gamma-ray bursts 
from neutron star collisions at cosmological distances (see, 
e.g., Katz & Canel 1996) if the energy in the fireball is suf- 
ficiently large, Er> ~ E~ l w E &na )> 10 51 <5f2/(47r) erg, and 
if the baryon loading of the fireball is sufficiently small, 
Mfb ^ 1CP 5 Mq (for a canonical energy Efb ~ 10 51 erg) so 
that the fireball can expand relativistically with a Lorentz 
factor Tfb = Ep a /{Mi\ } c 2 ) <; 100. Two questions arise 
from this suggestion. First, is the conversion of neutrino- 
antineutrino energy into electron-positron pairs enough 
efficient to provide the desired energy, and, second, how 
large is the baryon mass contained in the fireball created 
through vv annihilation? 

In the following we attempt to give answers to these 
two questions on grounds of our hydrodynamical collision 
models. Here we only report on the evaluation of the head- 
on collision Model H. We concentrate on this model for two 
reasons. On the one hand, the efficiency of vv annihila- 
tion increases strongly with the neutrino luminosities and 
with the mean energies of the emitted neutrinos. There- 
fore Model H with its larger neutrino emission appeared 
to us as the more interesting one. On the other hand, our 
simulations have demonstrated that even in the off-center 
collision the interaction of the two neutron stars is so dra- 
matic that a lot of matter is ejected perpendicularly to 
the orbital plane. Therefore, despite of the large angu- 
lar momentum in the system, the axis region is polluted 
with baryons (see Fig. |5|), and both the remnants of the 
head-on and off-center collisions adopt more or less spher- 



ical shapes after the dynamic interactions, with a central 
massive object being surrounded by a less dense, extended 
cloud of hot baryonic gas (Figs. |l] and ||). From this point 
of view Model O did not seem to offer better perspectives 
for the emergence of a relativistic fireball from a baryon- 
depleted region near the collided neutron stars. 

Figure gives a map of the energy deposition rate 
by vv annihilation into e + e~ pairs (averaged over the 
azimuthal angle around the z-axis according to Eq. ([[])) 
in the surroundings of the collision remnant in Model H 
at time t = 3 ms which is inside the double peak struc- 
ture of the neutrino luminosity of Fig. ^0|. One can see 
that the highest energy deposition rates of the order of 
3 x 10 30 erg cm~ 3 s _1 to 10 31 erg cm -3 s _1 occur immedi- 
ately outside the neutrinospheres but in layers with densi- 
ties still above and around 10 10 gcm~ 3 . At the displayed 
time, the integral energy deposition rate in matter with 
density below 10 11 gcm~ 3 (evaluated according to Eq. (^) 
is 3.4 x 10 52 ergs -1 . 

Since the neutrino luminosities and mean energies of 
the emitted neutrinos show significant variation during the 
computation time (see Figs. ^0[ pll and |22]) , we compute 
the time integral of the energy deposition rate, Eq. (^|), 
by employing the approximate treatment summarized in 
Eqs. (^)-(^). The phase between t « 1.5 ms and t « 3.5 ms 
has by far the highest neutrino luminosity and therefore 
yields the largest contribution to the time integral. For 
this reason we calculate the temporal average of the term 
L ann (t)/ T*(t) in Eq. (jj) by summing over N = 3 time 
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Fig. 26. Map of the local energy deposition rates 
(in erg cm" 3 s" 1 ) by vv annihilation into e + e" pairs in 
the vicinity of the merger for Model H at time t = 3 ms 
after the start of the simulation. The values are obtained 
as averages over the azimuthal angle around the z-axis, 
d measures the distance from the grid center in the x- 
y-plane. The corresponding solid contour lines are loga- 
rithmically spaced in steps of 0.5 dex, the grey shading 
emphasizes the levels with dark grey meaning high energy 
deposition rate. The dashed lines mark the (approximate) 
positions of the neutrinospheres of v e , v e , and v x (from 
outside inward), defined by the requirement that the op- 
tical depths in z-direction are r ZMi = 1. The dotted con- 
tours indicate levels of the azimuthally averaged density, 
also logarithmically spaced with intervals of 0.5 dex. The 
energy deposition rate was evaluated only in that region 
around the merged object, where the mass density is below 
10 11 gem" 3 . The integral value of the energy deposition 
rate at the displayed time is 3.4 x 10 52 ergs" 1 . 

points in this interval: t\ = 2.47 ms, ti = 3.01ms, and 
£3 = 3.59 ms. We obtain 

\ E (w) - 2 >< W- 57 MeV" 1 erg" 1 s . (10) 

i— 1 ^ V / / 

Because the three terms of the sum are different by less 
than a factor 2, we think that the splitting of the time 
integral of Eq. (0) which led to the approximate form of 
Eq. (|J) was justified. Taking the data for the individual 
neutrino luminosities L Vc (t), £ Pe (i), and E Vm (t) (Fig. f22| ) 
and the average neutrino energies {tv c (i)), (e P( ,(t)), and 
(e Vm (t)) (Fig. we further find 

plO ms 

/ F*(t)dt w 5 x lO^MeVerg^" 1 . (11) 
Jo 
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Fig. 27. Cumulative mass M(r > R) and annihilation 
energy E v „{r > R) outside of radius i?, as functions of R 
for Model H at time 3 ms (same time as in Fig. ^). The 
corresponding relativistic Lorentz factor T(R) = E uj} {r > 
R)/ [M(r > R)c 2 ] is also plotted. 

Multiplying the results of Eqs. ( [To| ) and (O) we end up 
with a total energy deposition of E ann w Kr° erg within 
our simulation interval of 10 ms for Model H (see also 
Fig. [27]). The corresponding average energy deposition 
rate by vv annihilation of ~ 10 52 ergs -1 is very large and 
means a conversion efficiency of vv energy to e + e" pairs 
of the order of 1%. Most of this energy is liberated during 
the 2 ms interval between t w 1.5 ms and t w 3.5 ms after 
the start of the simulation because of the enormous neu- 
trino luminosity shortly after the collision of the neutron 
stars. 

The energy of approximately 10 50 erg in e + e" pairs 
and photons is released nearly isotropically (Fig. ^6| ) and 
is only one order of magnitude below the canonical fire- 
ball energy En, » E 1 ~ 10 51 (5r2/(47r) erg. It may there- 
fore be sufficient to account for the shorter and weaker 
bursts whose energy is estimated to be typically more 
than a factor of 10 below the mean energy of the longer 
and more powerful bursts (Mao et al. 1994). However, 
most of the vv energy deposition happens very close to 
the neutrinospheres and thus in a high-density region 
(Fig. ^6|). Because of this, the baryon loading of the e + e"- 
photon fireball is a serious problem. This is obvious from 
Fig. ^ where we give the energy from ^P-annihilation 
and the baryonic mass, integrated from outside inward 
to the radial position R given on the abscissa. In the 
region where 10 50 erg are deposited, one has a baryon 
mass of M <; 5 x 10" 2 M Q which is about 5 orders of 
magnitude too large to allow for highly relativistic ex- 
pansion. The Lorentz factors which can be estimated as 
T(R) = E„v(r > R)J [M(r > R)c 2 ] are therefore around 
10" 3 instead of 100. For this reason, the huge and nearly 
isotropic baryon pollution of the vv energy deposition re- 
gion seems to rule out the possibility that neutrinos from 
colliding neutron stars produce gamma-ray bursts. 
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Fig. 28. The thick solid line gives the relativistic rota- 
tion parameter a(r) = J(r)c/(GAf 2 as (r)) for Model O as 
function of radius r. J(r) is the angular momentum (in z- 
direction) inside r, M gas (r) is the gas mass inside r which 
is plotted as thin solid line, given in units of Mq and scaled 
down by a factor 8. 

7. Summary and discussion 

We have reported results of three-dimensional Newto- 
nian hydrodynamical simulations of the collision along 
parabolic orbits of two identical, non-rotating neutron 
stars with a baryonic mass of about 1.6 Mq. The simu- 
lations were done with a Eulcrian PPM code, employing 
nested grids, using a physical nuclear EOS, and taking into 
account gravitational- wave emission and its back-reaction 
on the hydrodynamical flow as well as the emission of 
neutrinos from the heated neutron star matter. We have 
studied two different cases, each with varied resolution 
on the finest grid, a head-on collision and an off-center 
collision with a periastron distance of one stellar radius 
(ps 15 km). Our simulations allow us to draw the following 
main conclusions: 

1.) Dynamical evolution: During the head-on collision of 
the two neutron stars, the collision remnant is heated by 
recoil shocks and the initial kinetic energy is efficiently 
shock-dissipated to thermal energy within a few violent 
pulsations on a timescale of about 3-4 ms, after which the 
collision remnant has an essentially spherical mass distri- 
bution. The neutron star matter is transiently heated to 
peak temperatures close to 100 MeV. Average tempera- 
tures are around 40-50 MeV, corresponding to entropies 
between 2fc/nucleon and lOfc/nucleon. The violent crash 
of the stars into each other leads to the dynamical ejec- 
tion of about 0.5% of the system mass. In contrast, the 
off-center collision is much "milder" and the ejected mass 
is about a factor of 10 smaller, although the collision rem- 
nant retains a large fraction (about 30%) of the initial 
kinetic energy as rotational energy even after 10 ms at 



the end of our simulations. Only for extreme assump- 
tions about the nuclear EOS will the <; 3 Mq remnant 
of the head-on collision not form a black hole on a dy- 
namical timescale. We also believe that the remnant of 
the off-center collision will probably not be able to escape 
the collapse to a black hole within a few milliseconds. A 
mass of 3 Mq is distributed within a radius of 30 km and 
2.8 Mq are within 20 km (Fig. |28|) which is only about 
twice the Schwarzschild radius of a 3 Mq black hole. Ther- 
mal pressure can only insignificantly raise the maximum 
stable mass of neutron stars (Goussard et al. 1997) and 
rotation is able to increase the stable mass limit by only 
< 20% (Friedman et al. 1986). The interior - 2.8 Mq of 
the remnant rotate nearly uniformly and the relativistic 
rotation parameter of this mass is only a(r = 20 km) = 
[jc/{GM 2 )] 2Qkm « 0.3 (cf. Fig. H), which is smaller than 
that of the maximally rotating, maximum-mass models 
constructed by Friedman et al. (1986). For all but one ex- 
treme EOS tested by Friedman et al. (1986) such a config- 
uration is unstable. Therefore we conclude that it is very 
likely that gravitational instability will set in as soon as 
the two stars have merged into one massive body within 
t <; 3-4 ms. General relativistic effects will certainly influ- 
ence the transformation of the infall orbit of the neutron 
stars into a bound one; this needs to be studied with rel- 
ativistic simulations. 

2. ) Gravitational waves and neutrinos: The gravitational- 
wave signal will certainly depend on general relativistic ef- 
fects and our basically Newtonian models have only a lim- 
ited ability to make predictions of a possibly measurable 
pulse. Moreover, the duration of the gravitational-wave 
and neutrino emission from the hot collision remnant will 
depend on the timescale of the delay until black hole for- 
mation. Our simulations yield a maximum amplitude of 
the gravitational waveform that is ft. max ~ 2-4 x 10~ 23 for 
neutron star collisions happening at a distance of 1 Gpc. 
The off-center collision is the stronger gravitational-wave 
source due to the larger quadrupole moment of the ro- 
tating system and the longer duration of the emission 
which might last for 10-20 wave periods in the 1000- 
2000 Hz range. The gravitational- wave strain would be 
close to the lower sensitivity limit of the new generation of 
gravitational-wave interferometers which are currently un- 
der construction and will start operation within the next 
years. Of course, neutron star collisions are very rare and 
very short events and therefore the chance to catch a sig- 
nal is rather small. The head-on collision is the more pow- 
erful neutrino source of the two investigated cases and 
emits an energy about 3 times larger in neutrinos than in 
gravitational waves. The peak neutrino luminosity reaches 
4 x 10 54 ergs _1 and the total energy radiated in neutrinos 
within a few milliseconds is around 10 52 erg. 

3. ) Gamma-ray bursts: Because of the larger energy out- 
put in neutrinos and the very high neutrino luminosity 
(L„ > 10 54 ergs~ 1 ) as well as high mean energies of the 
emitted neutrinos ((e„) ^ 40 MeV), the head-on collision 
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provides more favorable conditions for producing gamma- 
ray bursts from e + e~-photon fireballs created by vv anni- 
hilation. We calculate a conversion efficiency of neutrino 
energy into e + e~-pairs of about 1% and find an integral 
value f or the energy deposited in the vicinity of the col - 
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lision fcmnant of 10 50 erg within only 10 ms , However, 

most of this energy is deposited in the immediate neigh- 
bourhood of the neutrinosphercs where the density is still 
higher than about 10 10 gcm~ 3 . Therefore the baryon load- 
ing of the e + e~-photon fireball is at least 5 orders of mag- 
nitude too high and instead of having values above 100 the 
relativistic Lorentz factor is estimated to be only around 
10 -3 . Dynamically ejected material together with a flow 
of baryonic matter driven by neutrino-energy transfer to 
the surface layers of the collision remnant are therefore a 
harmful poisonous combination which prevents relativistic 
expansion of the pair-plasma fireball even though the lat- 
ter seems to obtain an interesting amount of energy from 
w-annihilation. 

Strong shock heating of the neutron stars during their 
violent collision (or during their merging as suggested 
by post-Newtonian simulations, see Oohara & Nakamura 
1997) can indeed raise the neutrino luminosities signifi- 
cantly and can thus enhance the energy deposition by vv 
annihilation. The associated dynamical ejection of gas and 
the neutrino-driven wind caused by the intense neutrino 
fluxes, however, impede the emergence of gamma-rays 
from this scenario. Fireballs powered by neutrinos from 
collidir 
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: g iiuuuoii stars slioulu Uicicloie be n 
possib fe somcc'S of cosmolugical gaiiima-i &,y bmsts. Also, 

we do not see the formation of highly relativistic shocks 
during the collisions by which a fraction of 0.1-1% of the 
kinetic energy at impact (<; 10 53 erg) might be hydrody- 
namically focused into a small amount of mass (~ 10~ 6 to 
1O _5 M0). The expansion velocity of the collision shocks 
is at most a few tenths of the speed of light and the masses 
ejected in the presented models are more than 3 orders of 
magnitude too large to allow for ultrarelativistic motion. 
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